#ifndef __CGmom_CPP
#define __CGmom_CPP
#include "head.h"
#include "class_ma.cpp"
void cgmom(matrix_ &A,
	matrix_ &x,matrix_ &b,float err,COMPLEX x0)
//this function is based on Jin Jianming's"EM finite elemental method"
{
   int order=A.rows,k;
   matrix_ r(order,1),p(order,1),Ap(order,1),AHr(order,1),akp(order,1);
   Complex temp0;
   Complex ak,bk;
   double ru,ru0;
   double rss;
  
   
   ofstream fout;
   fout.open("iter_cg2_err.dat");
   fout.setf(ios::scientific);


   x=x0;
   //Ap=A*x;
   multi(Ap,A,x);
   //r=b-Ap;
   sub(r,b,Ap);
   //AHr=H(A)*r;
   Hmulti(AHr,A,r);
   //temp=H(AHr)*AHr;
   temp0=square(AHr);
   //p=AHr/temp[0];
   divid(p,AHr,temp0);
   //temp=H(b)*b;
   temp0=square(b);
   ru0=sqrt(abs(temp0));


   k=0;
	cout<<"\nstep of iteration="<<k;

   do
   {
      //Ap=A*p;
      multi(Ap,A,p);
      //temp=H(Ap)*Ap;
      temp0=square(Ap);
      ak=Complex(1./temp0);
      //x=x+ak*p;
      multi(akp,p,ak);
      add(x,x,akp);

      //r=r-ak*Ap;
      multi(akp,Ap,ak);
      sub(r,r,akp);

      //AHr=H(A)*r;
      Hmulti(AHr,A,r);
      //temp=H(AHr)*AHr;
      temp0=square(AHr);
      bk=1/temp0;
      //p=p+bk*AHr;
      multi(akp,AHr,bk);
      add(p,p,akp);
      //temp=H(r)*r;
      temp0=square(r);
      ru=sqrt(abs(temp0));
      rss=ru/ru0;
      //if(rss<err1)istodo=true;
      k=k+1;
      fout<<k<<","<<rss<<"\n";
		cout<<"\nstep of iteration="<<k<<";  error now="<<rss<<";  err="<<err;;
   }
   while(rss>err);
   fout.unsetf(ios::scientific);
   fout.close();

}
void cgmom(matrix_ &A,
	matrix_ &x,matrix_ &b,float err,COMPLEX x0,char *file)
//this function is based on Jing Jieming's"EM finite elemental method"
{
   int order=A.rows,k;
	matrix_ r(order,1),p(order,1),Ap(order,1),AHr(order,1);
   matrix_ temp(1,1);
   COMPLEX ak,bk;
   float ru,ru0;
   float rss;
   
   ofstream fout;
   fout.open(file);
   fout.setf(ios::scientific);


   x=x0;
   Ap=A*x;
   r=b-Ap;
   AHr=H(A)*r;
   temp=H(AHr)*AHr;
   p=AHr/temp[0];
   temp=H(b)*b;
   ru0=sqrt(abs(temp[0]));

   k=0;
	cout<<"\nstep of iteration="<<k;

   do
   {
      Ap=A*p;
      temp=H(Ap)*Ap;
      ak=1/temp[0];
      x=x+ak*p;
      r=r-ak*Ap;
      AHr=H(A)*r;
      temp=H(AHr)*AHr;
      bk=1/temp[0];
      p=p+bk*AHr;
      temp=H(r)*r;
      ru=sqrt(abs(temp[0]));
      rss=ru/ru0;
      //if(rss<err1)istodo=true;
      k=k+1;
      fout<<k<<","<<rss<<"\n";
		cout<<"\nstep of iteration="<<k<<";  error now="<<rss<<";  err="<<err;;
   }
   while(rss>err);
   fout.unsetf(ios::scientific);
   fout.close();

}



#endif